Steady-state and quench dependent relaxation of a quantum dot coupled to 

one-dimensional leads 
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We study the time evolution and steady-state of the charge current in a Single Impurity Anderson 
Model, using Matrix Product States techniques. A non equilibrium situation is imposed by applying 
a bias voltage across one-dimensional tight binding leads. Focusing on particle-hole symmetry, we 
extract current-voltage characteristics from universal low bias up to high bias regimes, where band 
effects start to play a dominant role. We discuss three quenches, which after strongly quench 
dependent transients yield the same steady-state current. Among these quenches we identify those 
favorable for extracting steady-state observables. The period of short time oscillations is shown to 
compare well to real-time renormalization group results for a simpler model of spinless fermions. 
We find indications that many body effects play an important role at high-bias-voltage and finite 
bandwidth of the metallic leads. The growth of entanglement entropy after a certain time-scale 
oc A -1 is the major limiting factor for calculating the time evolution. We show that the magnitude 
of the steady-state current positively correlates with entanglement entropy. The role of high energy 
states for the steady-state current is explored by considering a damping term in the time evolution. 

PACS numbers: 73.63.Kv, 73.23.-b, 72.10.Fk, 71.15.-m 



I. INTRODUCTION 



Over the past decade, experimental control over 
quantum systems has increased considerably. Possible 
realizations reach from model Hamiltoniansi>2 using 
ultra cold atoms in optical lattices to experimental 
setups of nanoscopic devices like molecular junctions, 
quantum wires or quantum dotsPS. Many of these sys- 
tems show remarkable properties, often due to reduced 
effective dimensionality and many body interactions. A 
prominent example is the Kondo effect^, which plays 
an essential role in transport across quantum dots. A 
theoretical understanding of transport in out of equilib- 
rium conditions is highly interesting for applications in 
nano- and molecular- electronics and even in biological 
systems. 

Electron-electron interactions render the theoretical 
description of non equilibrium dynamics one of the 
most challenging problems in today's condensed matter 
physics. However, with the advent of efficient numerical 
techniques to simulate one-dimensional (Id) quantum 



systems 
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many physical problems are well within 



grasp of theoretical physicists. 

In this work we obtain the steady-state charge current 
of a single interacting quantum dot under voltage bias, 
modeled by a single-impurity Anderson model (SIAM)ii. 
This model is commonly discussed in the wide-band limit 
approximation, tailored towards a universal, low-bias 
transport description. Here, we extend the discussion to 
the case of a finite (semi-circular) conduction band in the 
leads, which has not been explored specifically. A par- 
ticular realization could consist of two one-dimensional 
leads like nano-wire a 12 i 13 and a junction between them 
comprised of a magnetic impurity i.e. the quantum dot. 
We use generic one-dimensional tight binding leads with 



finite electronic bandwidth which mimic the electronic 
properties of for example carbon nano-tubes^. In such 
a device the electronic density of states (DOS) of the 
leads would have a bandwidth on the order of 15 e V 14 ' 15 
and effects arising from their specific structure are to 
be expected when using corresponding bias voltages 
which arc larger than those typically applied in current 
experiments with nanoscopic devices. 
The steady-state is obtained by combining Density 
Matrix Renormalization Group (DMRG)&±£ and Time 
Evolving Block Decimation (TEBD)^^ techniques, to 
perform real time evolution of the system after several 
different quenches. We focus on the particle-hole sym- 
metric point which shows the most pronounced many 
body effects^. We show that the same steady-state 
current is reached independent of the type of quench 
used and identify quenches which are superior to others 
when it comes to extracting steady-state data. We 
investigate quench induced oscillations in the transients 
and compare to real-time Renormalization Group re- 
sults. The method is well suited for reaching relevant 
time scales to study the steady-state current. We find 
that our approach is capable of yielding unbiased results 
valid in the thermodynamic limit. In the low bias region 
our results for the current-voltage characteristics agree 
with previous data (Heidrich-Meisner et al. Ref. Il7l ). 
We are able to extend earlier results^— to a larger 
parameter regime and discuss the interplay of finite lead 
bandwidth and electronic correlations. We find evidence 
for pronounced many body effects at high bias voltages 
in interplay with finite electronic bandwidths of the 
leads. Finally we discuss the role of high energy states 
for the steady-state current in low- and high bias voltage 
regimes. 

The text is organized as follows: In Sec. [H] we introduce 
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FIG. 1: (Color online) Illustration of the three quenches per- 
formed for the SIAM: i) QT I: quenching of both quantum 
dot-reservoir tunnelings t' L and t' R , ii) QT II: quenching the 
bias voltage Vb and iii) QT III: quenching the dot-lead tun- 
neling t' R to one lead. 



our model and describe in detail the different quenches 
to be performed. We present data for the transient 
response in Sec. Mil Results for the steady-state current 
are presented in Sec. |IV] where we also outline how to 
extract steady-state data from time evolved quantities. 
We analyze time scales of individual parameter regimes 
in Sec. |V] The role of high energy states in different bias 
regimes is discussed in Sec. I VII A detailed convergence 
analysis is presented in App. [5J 
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(see Fig.[I| where U parametrizes the on-site interaction 
strength on the quantum dot, t' a ,a £ {L,R} is the 
coupling strength between the quantum dot and the 
left and right lead. Lead a is characterized by intra 
lead hopping t and on-site potential e a . Particle- hole 
symmetry is enforced for all chosen parameters. When 
needed, the on-site energy of the quantum dot will be 
denoted by e/. 

We choose t = 1 and symmetric couplings 

all simula- 
4t of the 
widths of 

We will dis- 
We 



t' L = t' R = 0.3162 1 (Eq.flTd])) for 
tions, yielding a bandwidth of D = 
leads and an equilibrium Anderson 

Preservoir (0) ^ 0.1 £. 

play all energies in units of A (h, ks and e = 1 



restrict ourselves to the zero temperature case. Real 
time will be denoted by r. In App. [A] we show that 
within the simulation time r accessible, the fmitcness of 
the leads does not affect our results. 



II. SETUP 



In this section we define our notation for the SIAM 
which we use to model a quantum dot. We are interested 
in electron transport across the quantum dot after each 
of the several quenches to be described in detail in the 
following. We explain how we calculate the ground state 
using DMRG and the real time evolution using TEBD. 



A. Single-impurity Anderson model 



We consider a model for a quantum dot including 
charge as well as spin fluctuations described by the SIAM, 
consisting of an interacting site connected to a bath 
of non interacting electrons. We choose a setup where 
the quantum dot is located in the middle of a one- 
dimensional chain of tight-binding electrons. The dy- 



B. Quench preparation 

We are interested in the steady-state current^ of 
Eq. (fTaj) under a finite bias voltage Vb ■ Our strategy to 
obtain the steady-state is by quenching the Hamiltonian 
parameters xo = {U, t, t' al e a } at t = from some initial 
to their final values H(xo) — >• "H(x) and evolve an initial 
state \^o) with H(x.). \^>q) is chosen to be the ground 
state of the initial Hamiltonian H(xo) in the canonical 
ensemble at half filling with total spin projection S z = 0. 
We consider three different quench types (see Fig.[T]) 
which will be explained in detail below. Unless 
stated otherwise, we choose a system of L = 150 
sites with the quantum dot located at site 75. To 
drive the system out of equilibrium, a bias voltage 
Vb is applied by setting the respective on-site en- 
ergies of the leads in an anti-symmetric fashion to 
€l = —£r = -^p. For all quenches, the final parameters 
are x = {U,t = l,t' a = 0.3162 t,e L = ^f-e R = -^}, 
with variable U. The initial setup is quench type 
(QT)-dependent (sec Fig.H]): 
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1. QT I: Hybridization quench to both leads 
t' a = -> 0.3162 1 

For r < wc take x = {U,t,t' a = 0, e a = ±V B /2}, 
i.e. no quantum dot-to-leads coupling, but the bias 
voltage is already applied. We prepare the ground state 
of ri(x.o) at half- filling in the left and right lead and 
a single up-electron on the quantum dot. At r = 
the tunneling t' a is quenched to its finite value. Note 
that due to the splitting into three disconnected parts 
(t' a = 0), S z is not zero on the quantum dot and on the 
right lead initially. 



where N T = j- is the number of time slices, T is the 
total simulation time and St the length of a single 
time step. The operators 'K & and T-L act on even and 
odd bonds of the bipartite lattice respectively. Unless 
stated otherwise we use a TEBD matrix dimension of 
Xtebd = 2000 and a Trotter time step of St = 0.05 1 -1 . 
For additional details including studies of convergence 
in system size L and all auxiliary numerical parameters 
we refer the reader to App. |XJ 

The calculations carried out in this work set very high 
computational demands (« one million CPU hours) and 
were only possible due to a parallel code^r— which 
respects quantum number (TV, S z ) conservation. 



2. QT II: Quenching the bias voltage e a = — > ±Vb/2 

At r < 0, x = {U,t,t' a = 0.3261 t,e a = 0}. The 
system is prepared in the ground state \^o) at half filling 
with overall S z = zero. At r = the bias voltage is 
quenched to its desired value. As compared to QT I, this 
setup has the advantage that no subsystems with finite 
values of S z exist in the ground state. Furthermore, 
correlations between the three regions are already 
present in the ground state. Note that the initial state 
is much more complicated than for QT I. This typ e of 
quench has also been used by the authors of Ref. \24 . 



3. QT III: Quenching the hybridization t' R — — > 0.3162 1 
to the right lead 

The initial parameters are chosen Xo = {U,t,t' L = 
0.3261 £, = 0, e a = ±Vg/2}, and the system is again 
solved for the ground state | \&o ) at half filling. At r = 0, 
we quench t' R , = — >• 0.3162 1 and evolve \^q) with the 
quenched Hamiltonian. An advantage over quenching 
both quantum dot-reservoir tunnclings (QT I) is that 
the two disconnected parts in this approach initially 
have S z = 0. 



C. Methods 

To prepare the system in the ground state of the ini- 
tial Hamiltonian, we employ the DMRG^ algorithm in 
its two and single site formulation. Our implementa- 
tion exploits conservation of spin projection (S z ) and 
charge (N), which is crucial for obtaining high precision 
data. Time evolution is done using the TEBD 8 algo- 
rithm, within a second order Suzuki- Trotter decomposi- 
tion of the propagator 

e-^ T = U-™ 5t )^ = (eT^e^'eT"')"' + O (St 3 ) 



III. TRANSIENT RESPONSE 

In this section we present results for the transient cur- 
rent response of the three quenches. We discuss individ- 
ual bias regimes and identify oscillations in the time evo- 
lution of the current which are reminiscent of results for 
an interacting resonant level model of spinlcss fcrmions. 
We show that QT II leads to much larger initial oscilla- 
tions than the other two QTs. 



A. Low, medium and high bias regime 

In our simulations we identify three regimes of bias 
voltage Vb with qualitatively different behavior. Within 
each regime, the general features of the time evolution 
of the current are only moderately dependent on inter- 
action strength. For that reason, we first discuss results 
for U/A = 12 only (sec Fig.©. 

For low bias voltages (Vb/A £ (0,18)), a steady-state 
current plateau is reached within r ps A -1 . 
In a region of medium bias voltages (Vb/A £ (18,28)) 
we observe a fast increase in current over a timcscale of 
r rs 0.3 A -1 followed by a rather slow decay which, for 
some model parameters, is too slow to reach a steady- 
state plateau within accessible simulation times (see be- 
low). 

In the high bias region (Vb/A £ (28,40)) the time evo- 
lution of the current shows a sharp peak followed by fast 
decrease of the current into a steady-state plateau within 
rw A- 1 . 

Our data indicates that within a simulation time of 
r = 3 A -1 approximately one particle is transferred from 
the left reservoir to the right one. As discussed in detail 
in Sec. IIV| all three QT eventually approach the same 
steady-state, although in quite different manner. QT II 
for example leads to the largest transient current spike, 
which is one reason for the lower accuracy in determin- 
ing the steady-state for this quench. We also find that 
quenching the hybridization(s) (QT I or III) yields much 
> cleaner steady-state plateaus as compared to quenching 
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FIG. 2: (Color online) Time dependence of the current (Eq. (|A1[) ) at U/A — 12 for the three different QTs and for different 
bias voltages. The curves are plotted as solid lines up to the last reliable point in the TEBD calculation (see text). Larger 
times are plotted as dash-dotted lines. Solid horizontal lines are fits to extract the steady-state currents. The time domain for 
these fits starts at r ~ A -1 and ends at a point which is identified as the last reliable data point (symbols, see text). Dashed 
horizontal lines indicate the uncertainty. The insets in the mid row show respective zooms onto short time regions, which are 
not visible in the main part of the figure. 



the bias voltage (QT II), which leads to more pronounced 
oscillations in these plateaus. 



B. Characteristic oscillations of the current 

The time evolution of the current exhibits oscillations 
which are more or less pronounced depending on the type 
of quench. These oscillations become more explicit with 
increasing interaction strength (not shown) . Their period 
is of the order of 0.5 A -1 for low bias voltages, increasing 
to about 0.3 A^ 1 for higher bias voltages. These oscil- 
lations compare nicely to results from real-time renor- 
malizationgroup for the interacting resonant level model 
(see Ref. [28|, equation 107), which predicts a sinusoidal 



behavior (oc sin(^)) with a period of 



tc(U,V) = 



V B + U 



In Fig. |31 we plot Tc(U,V) as a function of interac- 
tion strength and find remarkable agreement with rtRG- 
results at higher bias voltages. The period was extracted 
from the TEBD time evolution data in three ways: i) 
by fitting a sine function, ii) by identifying the domi- 
nant Fourier amplitudes and iii) by identifying local max- 
ima. These data were combined and their individual un- 
certainty taken into account. Error bars in Fig. [3] (not 
shown) would be sharply growing below Vb = 25 A. In 
the lower bias regime our data is not significant for a re- 
liable extraction of the period. 
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FIG. 3: (Color online) Period of the sinusoidal oscillations of 
the current in QT II for various values of interaction strength 
U/A = 0,4,8, 12 and 20 (symbols). Solid lines indicate the 
predicted form for the interacting resonant level model2£. 
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FIG. 4: (Color online) Current-voltage characteristics of the 
quantum dot. The steady-state currents shown are obtained 
by a fit of the expectation value of the current operator within 
the steady-state plateau. Regions where only an upper bound 
for the steady-state current could be obtained are indicated 
by pedestals (see text). 



IV. STEADY-STATE CURRENT 

In this section we present the current-voltage charac- 
teristics of quantum dot. We outline a scheme to extract 
the steady-state current and investigate the dependence 
on the type of quench used. The current- voltage charac- 
teristics in the low bias region is compared to existing 
data obtained with other methods. Furthermore we 



present a detailed comparison between an interacting 
and a non interacting quantum dot for finite as well as 
infinite lead bandwidth. 



A. Extracting the steady-state current 

We identify the steady-state current as the mean value 
of the time-dependent current taken over a suitable time 
domain [ts,te] over which the current shows an almost 
constant behavior (apart from small oscillations). t$ 
typically depends on the model parameters and was 
chosen by hand, and te is taken to be the largest time 
for which simulations yield reliable results (see Fig. [2]). 
Beyond te the current becomes numerically unreliable, 
resulting in an artificially decaying current (see App. [S] 
for discussion). We find that in most of the parameter 
regions the transients have decayed at t$ ~ A -1 . On 
the other hand, the end point of the plateau strongly 
depends on the parameter region under consideration. 
We define it by two distinct measures. One is the time 
for which the truncated weight e (see Eq. (|A3|) ) 
reaches a threshold of e c = 3 • 10~ 6 at any bond (marked 



by + in Fig. [2}. The second definition (r|?' marked by 
o in Fig. [2]) is given by the time for which two different 
definitions of the current, namely the expectation value 
of the current operator (Eq. (|A1[) ) and the time deriva- 
tive of the particle number (Eq. (|A2|) ). deviate by more 
than 7 • 10 -4 , the latter being more susceptible to accu- 
mulation of errors. Both times are in good agreement 
with each other and can be combined into an effective 



simulation time te 



■ I (i) 
mm(T E 



.(21 



) + «| 



.(2) 



(marked by triangles in Fig. [2]). We choose a value of 
a = 0.1. Results do not depend on this particular choice. 
The steady-state plateaus obtained in this way usually 
show oscillations and/or small, parameter- and quench- 
dependent drifts. We quantify the quality of convergence 
within the plateau region [ts, te] by the slope of a linear 
fit to the current. A large slope indicates that it is 
not possible to reach the steady-state within the given 
simulation time te, i.e. the physical relaxation time is 
too long or the reached simulation time is too short. 
This is further discussed in Sec. [Vj For these parameter 
values we can only extract an upper bound for the 
steady-state current, given by the current at the last 
reliable simulation time. This is justified because we 
find the current to always decrease as a function of 
time (apart from small oscillations). We consider the 
current to be converged when the relative slope is below 
a threshold of ~ 5 ■ 10~ 2 A. Each curve in addition was 
inspected by hand for convergence. When we consider 
the steady-state current converged, we estimate its 
error as three times the standard deviation taken over 
the data points in the fitting interval [t$ , te] (plotted 
as dashed lines in Fig. [2] and Fig-HJ) - This coincides 
most of the time with the maximal deviation of the 
time-dependent current from its mean value. 
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As an important test, we obtained the current-voltage 
characteristics for the non interacting case and compared 
it to analytical results^ (see Fig.[4j) , finding excellent 
agreement. Another indication for the reliability of the 
scheme outlined above is that all three types of quenches 
investigated yield the same steady-state current within 
the uncertainty. As noted in App. |A] the position of 
the plateau is also stable with respect to variations of 
technical parameters of the simulation. The quality of 
the steady-state plateau, however, depends strongly on 
the values of interaction and bias voltage and may be 
obscured by initial oscillations or shortened at the end 
by the truncated weight breakdown. 
The behavior of the spin-current strongly depends on 
the quench type and is even identical to zero for QT 
II. In this respect, the steady-state charge current does 
not depend on the properties of the spin current, since 
all three quenches yield the same steady-state for the 
charge current. 

From our calculations, we find QT I and QT III to yield 
more reliable data for the extraction of the steady-state 
current than QT II. Reasons for this behavior are i) the 
much more pronounced oscillations in the data of QT II 
which enlarge the statistical uncertainty of steady-state 
values and ii) the much higher transient spike in QT 
II which causes a higher entanglement and shortens 
te- In the following, we will present steady-state data 
extracted from QT I and QT III. 



B. Current-voltage characteristics 

The current-voltage characteristics of the quantum dot 
for interaction strengths of U/A — 0,4,8, 12 and 20 arc 
shown in Fig. [4] We plot data as obtained from QTs I 
and III (other QTs would give the same results but with 
larger error bars, as discussed in Sec. IIV A|) . Results for 
the non interacting case agree with analytic results for an 
infinite system 42 . In some regions only an upper bound 
for the steady-state current can be obtained. This region 
does not lie on the extreme end of the parameter space. 
It shows non trivial dependence on U and Vb, which is 
discussed in detail in Sec. [Vj The current- voltage charac- 
teristics has an approximately semi-circular shape, with 
decreasing maximum as a function of interaction strength 
U. At small bias Vb, the current is linear in Vb and 
agrees with the linear response result j lin = 2G Vb (see 
also Fig. [5]). At higher bias, it departs from the linear 
response result. With increasing U , this departure oc- 
curs already at smaller bias Vb , which can be attributed 
to an exponential thinning of the Kondo resonance with 
increasing U . 

In intermediate bias regions, we observe a flattening in 
the current-voltage curve. The maximum steady-state 
current is obtained in a bias regime from Vg ~ 15 A to 
Vb ~ 19 A, but is non monotonic as a function of U. 
Increasing the interaction from U = to U = 12 A shifts 



the maximum position to higher bias, whereas for higher 
U the maximum position again moves to lower bias volt- 
ages. 

The decrease of the steady-state current for high bias- 
voltages can be attributed to the diminishing overlap 
of the DOS of the two reservoirs. Both have a semi- 
circular DOS with a bandwidth of D = 40 A. In the 
wide-band limit, the curves behave similarly inside the 
low bias regime but should saturate as a function of Vb 
for higher bias voltages (see Fig. [5]). 
We discuss three simple limits. The TEBD results for the 
current respect the linear response (jii n ) for very low bias 
voltages which gives the conductance quantum Go- Fur- 
thermore they respect the high bias voltage band cutoff 
where the current has to go to zero (here at Vb = 40 A) 
due to diminishing overlap of the DOS of the reservoirs. 
The third limit is the non interacting case (non-trivial 
for the used numerical method), where we obtain perfect 
agreement with analytical results for the thermodynamic 
limit. 



C. Comparison to previous results 

In the low bias region, results from other techniques 
are available for the SIAM out of equilibrium. In the 
following we discuss our results for various values of 
interaction strength U together with data previously 
obtained (see Fig. [5]) by diagrammatic QMCi£, fourth 
order Keldysh perturbation theory^, time-dependent 
DMRGil, TEBD for temperature T = (this work), 
non equilibrium FRG^, non equilibrium Cluster Per- 
turbation Theory 2 ^, the non equilibrium Variational 
Cluster Approac h 21 ' 29 , imaginary time QMC2 2 -, iterative 
summation of real-time path integrals 2 ^ and the linear 
response result for the Kondo regime ju n = 2GqVb- All 
methods work at or close to zero temperature. Some 
of the methods use a wide-band limit and others a 
semi-circular reservoir DOS which (for equal A) become 
comparable in the shown low bias region (see Fig.|B] for 
a comparison). The newly obtained TEBD results agree 
very well with the unbiased dQMCi£ and quasi-exact 
tDMRGil data. An earlier comparison including more 
details but fewer techniques is available in Ref. |3Cj. 



D. Comparison to a non interacting device: 
Identifying correlation effects from the steady-state 
charge current 

To gain further understanding of the role of correla- 
tions, we compare the steady-state current of the inter- 
acting quantum dot (U, on-site potential e/ = —U/2) 
to the one of a corresponding non interacting (resonant 
level) device with U = and on-site potential e/ = — ¥ 
(sec Fig. [5]). Data for the resonant level device are ob- 
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FIG. 5: (Color online) Comparison of the current-voltage 
characteristics of the SIAM obtained with different methods 
in the low-bias regime. Some of the methods use a wide- 
band limit and others a semi-circular reservoir DOS which (for 
equal A) become comparable in the low bias region shown. 
The methods are: 1) diagrammatic QMC for T = in the 
wide-band limit (dQMC)^, 2) fourth order Keldysh pertur- 
bation theory for T = in the wide-band limit (PT4)±£ 3) 
time-dependent DMRG for T = using a semi-circular DOS 
(tDMRG)i£, 4) TEBD for T = using a semi-circular DOS 
(TEBD, this work), 5) non equilibrium FRG for T = us- 
ing a wide-band limit (FRG)— ^, 6) non equilibrium Cluster 
Perturbation Theory for T = using a semi-circular DOS 
(nCPT 11 )— , 7) non equilibrium Variational Cluster Approach 
for T = using a semi- circular DOS (nCPT^, 8) imaginary 
time QMC for T = 0.2 A in the wide-band limit (cQMC)^ 9) 
iterative summation of real-time path integrals for T = 0.2 A 
in the wide-band limit (ISPI)-^ and 10) the linear response 
result for the Kondo regime j\ in — 2GqVb (lm. resp.). 



tained analytically . 

From the plots in Fig. [HI one can sec clear differences in 
the low bias region between the non interacting and inter- 
acting device for all interaction strengths, which can be 
attributed to the presence of the low energy Kondo res- 
onance in the interacting case. For low bias, the Kondo 
resonance fixes the linear response current to a U inde- 
pendent constant and causes a higher current than for a 
non interacting quantum dot at the same on-site poten- 
tial. Furthermore the curvature of the current-voltage 
characteristics in the low bias region is negative in the 
interacting case as compared to positive in the non inter- 
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FIG. 6: (Color online) Comparison of the current-voltage 
characteristics of a non interacting, resonant level device with 
on-site potential €f — — y (solid lines) with the TEBD data 
for the interacting quantum dot (symbols) . Both devices have 
the same specifications with only the interaction U missing in 
the first case. The comparison is done for four values of inter- 
action strengths resp. on-site potentials: ^ = {4, 8, 12, 20} 
resp. -£ — {—2,-4,-6,-10} (blue/circles, green/triangles, 
red/stars, cyan/squares respectively). In addition we show 
the (7 = result (black/no symbols). The dash-dotted lines 
indicate data for a non interacting device in the wide-band 
limit. 

acting system. For larger values of U = 12 A and 20 A, 
the negative curvature turns into a positive one in the 
low bias region. 

For low values of interaction strength (see data for U — 
4 A) we observe deviations in both, the low and high 
bias region. For the latter, this hints at possible many 
body effects which may also be important in the high 
bias regime. Data in the medium bias region are almost 
indistinguishable from the non interacting case. For high 
values of interaction strength the picture changes and 
many body effects are present in the whole bias regime. 
Summing up, we find that effects of interaction are most 
pronounced in the low and also in the high bias regime, 
where a larger current is obtained than in the non in- 
teracting device. Because of the small remaining overlap 
of the DOS of the leads this larger current may be due 
to some low energy spectral weight in the interacting de- 
vice, consistent with low energy excitations observed in 
Ref. Hi] using a non equilibrium Variational Cluster Ap- 
proach calculation. 



V. DISCUSSION OF TIME SCALES 

We start out by discussing Kondo correlations. We 
then note that depending on bias voltage and interaction 
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strength, the steady-state charge current cannot always 
be reached within the reachable simulation time te (see 
Sec. IIV A[) . which can be attributed to i) weak spots 
of the method (i.e. small te) and/or ii) long physi- 
cal relaxation times. To obtain insight into physical 
mechanisms as well as the parameter dependence of the 
performance of TEBD, also relevant for future studies, it 
is desirable to disentangle these two effects. We identify 
parameter regimes with long physical time scales to 
be at U + Vg > D, which is a different parameter 
dependence than for the weak spots of the method, 
which are associated with large values of the current 
magnitude, as shown in App. [Bj 



A. Finite simulation size/time and Kondo 
correlations 

At the particle-hole symmetric point of the SIAM, 
Kondo correlations arc especially pronounced. In equi- 
librium they introduce a characteristic energy scale, the 
Kondo temperature which translates into a length 
scale of the Kondo singlet 6f which is given by Bcthc 
Ansata^i 



6f oc 



VF 

^bTk 



oc 2t 



AU 



(2) 



Due to the exponential dependence on interaction 
strength, these spin correlations can not fully develop on 
a finite size syste m 32 i 33 , already for moderate interaction 
strength. For the parameters used in this work the 
equilibrium Kondo correlations have a spatial extent 
(screening cloud) of approximately 6f ~ 50 sites for U = 
4 A, 6r ~ 200 sites for U = 8 A, £ K « 900 sites for U = 
12 A and £ K « 16000 sites for U = 20 A (see Eq. ©). 
Comparing QT I (disconnected before quench, no 
Kondo correlations) and QT II (leads connected, Kondo 
correlations present in initial state), we find the same 
charge current after a short transient regime. We 
conclude that although finite size systems are not able 
to capture the full equilibrium Kondo singlet 3 ^, the 
steady-state transport in the charge channel is not 
noticeably affected. Since Vb > Tk in all simulations, 
we also do not expect issues due to the Kondo time scale 
during the non equilibrium time evolution. 



B. Time scales in the high bias regime 

We find that relaxation times in the model under 
discussion are strongly parameter (U, Vb) dependent. 
These relaxation time scales are estimated by the slope of 
a linear fit to the plateau region [ts,te] (see Sec. IIV A[) . 
In particular, we identify three regions (see Fig. [7] (top)): 
region I is characterized by short physical relaxation 
times and region II exhibits longer relaxation times. 
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FIG. 7: (Color online) (top) Parameter regions in the U — Vb 
of short (I) and long (II) physical relaxation scales as well as 
a regime of more complex behavior III. Data from the TEBD 
calculation is indicated with black and grey markers. For 
region II, pedestals are shown in Fig.|3] (bottom) Single par- 
ticle DOS and single particle dot level in a Hubbard-I type 
picture at U = 20 A, for (a) V B = 6 A, (b) Vb = 20 A and 
(c) V B = 36 A. The electronic DOS of the left (right) lead 
is shown in red (blue) and their overlap in brown. The sin- 
gle particle level of the quantum dot is indicated at — ^ in 
magenta. 



Region II overlaps with the regime in which TEBD 
restricts us to small final simulation times te (high 
steady-state current regime, see App. |B]for discussion). 
In region II we did not obtain a converged steady-state 
current. In region III, the current is small and the 
maximum reachable simulation time (see App. [B]) was 
large enough to determine the steady-state current. 
We proceed by providing an intuitive single-particle 
picture of the transition from region I to II in a Hubbard 
I type description (Fig. [7] (bottom)). Then the leads 
(assuming infinite reservoirs) arc described by semi- 
circular bands of bandwidth D, asymmetrically shifted 
against each other with increasing bias voltage Vb- The 
quantum dot consists of a single (non interacting) level, 
located at the single particle energy — ~. We find that 
the transition occurs when this single particle level of 
the quantum dot leaves the overlap region of both lead 
DOS (blue line, U trans ~ D — Vb)- We conclude that 
the existence of an appreciable spectral weight in the 
overlap region of the lead DOS leads to faster relaxation. 



9 



VI. ROLE OF HIGH ENERGY STATES 

To study the role of high energy states during the time 
evolution we add a damping term to the propagator 



1.5 



U(r) 



_ -mT(i-tr) 



(3) 



which gradually reduces the contribution of high energy 
states. 

In Fig.0 the effects of damping of high energy modes 
on the current are visualized. We show results for very 
low bias voltage (Vb = 2 A) as well as high bias voltage 
(Vb = 32 A). The different influence of over-damping 
(dashed lines) on low bias setups in contrast to high bias 
setups yields insight into the role of high energy states in 
the two respective cases. In low bias settings, strong over- 
damping (here T = 10 A) leads to lower current while in 
the high bias case it leads to higher current with respect 
to the true one. This indicates a qualitatively different 
role of high energy states for these two settings. 
This result can be made plausible by a simple argument. 
In the case of small bias voltage (Vb << t), the dominant 
energies should be the kinetic ones and neglecting high 
energy states amounts to eliminating those with highest 
kinetic energy. Such states contribute much to the cur- 
rent and neglecting them leads to a lower total current. 
On the other hand for very high bias voltage (Vb »t), 
potential energy is expected to dominate. High energy 
states are then those with a lot of particles in the high 
bias reservoir. Eliminating them reduces the available 
state space for hopping of particles back to the side of 
high potential. Therefore the current is increased due to 
less back flow. 

From a technical point of view, such an approach may re- 
duce entanglement growth (the limiting quantity in real 
time evolution using matrix product states), thus reduc- 
ing the required matrix dimensions of the MPS. Using 
such an Ansatz however suffers from two drawbacks, i) 
On the one hand, we have just seen that high energy 
states can be important for the steady-state current, and 
on the other hand, estimating a priori a suitable mag- 
nitude of the damping T is not straightforward, since it 
should in principle be dynamically adjusted during the 
time evolution taking into account energies and trun- 
cated weight. Due to these reasons we refrain from using 
such an approach in general. However we show in Fig. [8] 
that by choosing a phenomenologically good value for the 
damping (T = A), one can indeed somewhat prolong the 
stable time evolution. 



VII. CONCLUSIONS 

We studied the single-impurity Anderson model out 
of equilibrium beyond the linear response regime by 
means of Density Matrix Rcnormalization Group. Real 
time evolution was performed making use of the Time 
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FIG. 8: (Color online) Effects of damping T of high energy 
modes on the time evolution of the current (U = 0, Xtebd = 
500, QT I). Data shown are obtained for very low bias volt- 
age (Vb = 2 A, group of grey curves in lower part of figure) 
and high bias voltage (Vb = 32 A, group of orange curves in 
upper part of figure). For each bias voltage we compare data 
obtained by a standard (r = 0) time evolution (full lines), 
data using an (empirically) optimally damped time evolution 
(r = A, dash-dotted) as well as for an over-damped evolution 
(r = 10 A, dashed). 



Evolving Block Decimation algorithm which allows to 
access relevant time scales to reach the steady-state. 
Within this framework we investigated three different 
quenches: i) quenching the hybridization with already 
applied bias voltage, ii) quenching the bias voltage and 
iii) quenching the hybridization at one side only. 
Calculated current-voltage characteristics agree very 
well with established results which are available in the 
low bias region. We find that the period of characteristic 
oscillations in the time evolution of the charge current 
is already very well described by rcnormalization group 
results for a different model, the interacting resonant 
level model of spinless fcrmions. After an initial 
transient regime, where on the order of one particle is 
transferred through the quantum dot, the steady-state 
current agrees among the three quenches investigated. 
For the identification of steady-state plateaus in time 
dependent quantities the type of quench is however 
very important. We show that quenching the lead-dot 
tunnclings is the most suitable one, contrary to expec- 
tations whereas quenching the bias voltage results in 
large initial oscillations of the current. We furthermore 
show that limitations of the method like its inherent 
finite size do not pose a problem for simulations of 
the setup discussed here within accessible times. Our 
findings indicate that the steady-state charge current 
is not influenced by finite size effects, hinting that 
incompletely developed Kondo correlations in the spin 
channel do not influence charge transport noticeably. 
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We find that a large entanglement entropy correlates 
positively with a large steady-state current amplitude. 
By studying a damped time evolution we find that high 
energy states have very different significance in the low 
and high bias regime respectively. 

Besides reproducing the universal low bias physics we 
open up new perspectives for devices in which a large 
bias voltage is combined with a finite electronic DOS 
of the leads, like nano-tubcs. For such devices we 
predict that effects of electron-electron interactions are 
important even at high bias voltages. 
Interesting extensions within the presented approach 
may be the application of a gate voltage to study 
stability diagrams, evaluation of spin correlations which 
could hint on Kondo correlations, to study effects of 
asymmetric couplings, the interplay of bias and magnetic 
fields as well as to investigate correlated leads. On the 
technical side it would be interesting to evaluate whether 
more gently ramped quenches over a finite time interval 
further decrease oscillations or even entanglement and 
further improve the extraction of steady-state data. 



Robinson bound^i, up to exponentially suppressed parts, 
and in our case is v w 2t, which can also be clearly 
seen in the time evolution of local charge expectation 
values 3 ^. The perturbation will hit the left and right 
end of the chain and return back to the quantum dot 
after a time of about r w 2(L/2)/(2£) = L/(2t), i.e. 
t/A' 1 fn LA/(2t) w 7.5 for L = 150. This is far beyond 
the times te (see Sec. IV A) up to which we calculate 
the steady-state current, which is therefore not affected. 
This conclusion is confirmed by the convergence of the 
current with respect to system size L, discussed below. 
The measured current may however be affected by other 
possible errors within our approach: i) the procedure 
to measure it, ii) the Trotter error, and iii) the limited 
matrix dimension xtebd (i.e. truncated weight). 
In the following we will show that the major uncertainty 
arises from the limited matrix dimension xtebd, while 
other source are negligible. The definition of the time 
intervals from which the steady-state current is evalu- 
ated (Sec. IV A) is also relevant. A similar conclusion 
has been drawn before in the framework of adaptive 
tDMRG^i. 
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Appendix A: Method setup, convergence analysis 
and preliminary considerations 

Here we present some preliminary considerations 
concerning the convergence and quality of our data. 
Uncertainties arise from the approximations made within 
the method and from numerical precision. 
In addition, our setup contains leads of finite size, with 
two effects in principle. First, this finite size affects the 
ground state at time zero. We will show below that the 
effect on the current is negligible; it converges already 
at much smaller lead size than used here. 
Secondly, the finiteness of the reservoirs means that 
no energy or particle dissipation occurs and eventually 
the system will show oscillatory behavior. We note 
in passing that during our simulation time only ap- 
proximately one particle traverses the quantum dot. 
The earliest time at which the current can be affected 
arises from a perturbation which propagates after the 
quench to the end of a lead and back to the quantum 
dot. The velocity of this signal is limited by the Lieb 



1. Obtaining the current 

Within the TEBD time evolution the steady-state cur- 
rent may be obtained via the expectation value of the 
current operator at each time step 36 ' 37 

3ij i T ) = 1 Uj ( a L a ja - a ia aj CT ) , 

where i and j denote adjacent sites and ai a and a ia are 
annihilation and creation operators for fermions on-site 
i with spin a which depend at time r and is taken to 
be real. To obtain the current though the quantum dot, 
a symmetrized version of the inflow and outflow is used 




where c^ nda , c;^ , denote operators on the last site of the 
left reservoir (number 74 in Fig.[l} and Cg^,cj^ denote 
operators on the first site of the right reservoir (number 
76 in Fig.p. 

Another way of computing the current is by calculating 
the time derivative of the total particle number to the 
left of the site under consideration 

3u+i(t) = ^ / ^2 S hm ^ T ) J ■ 

\ni— 1 cr / 
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FIG. 9: (Color online) Convergence of the current with respect to several auxiliary numerical parameters, (left) Full lines denote 
results obtained evaluating the expectation value of the current operator Eq. (| A If) . while dashed lines indicate data obtained 
by evaluating the time derivative of the expectation value of the particle number Eq. (|A2[) (U = 12 A, L = 150, Xtebd, 
QT I), (center) Matrix sizes xtebd = 250 (dotted), 500 (dash-dotted), 2000 (dashed) and 4000 (dash-dotted) are presented 
(U/A = 20, L = 150, QT I), (right) We show system sizes L = 20,40,60,80,100,120 and 150 (dotted, dash-dotted, dashed, 
dash-dash-dot-dotted, long-dash-short-dashed, dash-gap-dashed and solid) at U = 0, L = 150, Xtebd = 2000 for QT I. 



Again a symmetric combination of the dot's in- and out- 
going current yields the current under consideration 




The current through the dot may be evaluated at each 
TEBD time step using Eq. (|A1[) or by computing a finite 
difference approximation to the differential Eq. (|A2[) 
every two successive time steps. Besides the expected 
additional source of error by evaluating the time deriva- 
tive numerically, this method is expected to perform less 
well due to the influence of all sites in the system on the 
result for the current. Since the occupation number for 
each site has its own limited accuracy. A comparison of 
the current evaluated by means of Eq. (|Alj) and Eq. (|A2|) 
for various values of interaction strength U and applied 
bias voltage Vb as well as all QTs (I, II, III) shows good 
agreement in the beginning of the time evolution (see 
Fig-IU (left)). Due to an accumulation of errors in the 
particle number expectation values of the individual 

(2) 

sites, the results start to deviate at some time t e . We 
do not use results beyond t e (see the discussion in 
Sec. IIV Al Numerical values of all steady-state currents 
will be obtained using the current operator (Eq. (|A1[1 ) 
which yields a much more stable estimator. 



2. Finite size effects: L 

In this section we discuss the dependence of the 
results for the current on the length of the system L. 



We quench both dot-lead tunnelings (i.e. QT I). The 
qualitative behavior for the other QTs (II and III) is 
virtually identical. Results for the steady-state current 
for system sizes of L = 20, 40, 60, 80, 100, 120 and 150 
sites are shown in Fig. [5] (right) in the non interacting 
case. We find that the final results for the steady-state 
current agree with the analytically available results for 
an infinite system in all cases. This ensures a reliable 
determination of steady-state properties even on finite 
size systems. As mentioned before the system size 
limits the maximum simulation time due to signals 
back-propagating from the borders. In the following 
all calculations will be performed for a system size of 
L = 150 to provide a nice long plateau (maximum 
simulation time) in the steady-state current. 



3. Trotter error: St 



Usually the contribution to the total error arising due 
to the Trotter approximation is negligible with respect 
to other approximations. Since i) The Trotter error 
grows only linearly with simulation tim o 24 ' 38 , and ii) it 
can be controlled by choosing sufficiently small St. Wc 
investigated the influence of the Trotter decomposition 
on the current. Results for Sr/t^ 1 = {0.01,0.05,0.1} 
were found to agree to within 5 • 10~ 5 . We do not plot the 
results because they all lie on top of each other. A good 
value for the time step was found to be Jr/i" 1 = 0.05 
which was used in the main section of the paper. 
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4. MPS matrix dimension: \ 

The quality of the TEBD results is predominantly de- 
termined by the maximum matrix size Xtebd used. A 
bigger xtebd leads to fewer discarded states (i.e. less 
truncated weight of the reduced density matrices) during 
the truncation and therefore to a better approximation. 
The truncated weight is defined by2& 



size of xtebd = 2000, a trotter step of <5t/A _1 = 0.5 
and evolving for 1000 time steps which yielded a final 
simu lation time of T/A^ 1 = 5. Requiring a maxim um 



7=1 



(A3) 



where denote the eigenvalues of the reduced density 
matrices. This quantity is zero if no truncation is done. 
The computational cost of the TEBD algorithm scales 
essentially like^ 

cost cx L(d xtebd) 3 , 

where L is the length of the chain and d = 4 the size 
of the local fermionic Hilbert space. Therefore it is 
essential to keep xtebd as low as possible. During 
the simulations we noticed that at a certain time (long 
before signals propagating back from the ends of the 
chain would reach the quantum dot) the truncated 
weight starts to grow quickly and the results become 
unstable, causing a decaying current. The effects of 
enlarging xtebd are shown in Fig. [5] (center). As the 
data indicates, the effect of increasing xtebd is to make 
larger simulation times accessible, before the simulation 
breaks down due to accumulation of truncated weight. 
Remarkably, no spurious quasi steady-state is entered 
when Xtebd is relatively small. The overall shape of the 
current appears to be unaffected by enlarging xtebd, 
making reliable predictions for xtebd = 2000 possible. 
Therefore in the main part of the paper we always used 
Xtebd = 2000 as a good compromise between run time 
and accuracy. 



5. Setup 

Based on the above considerations all data in the 
main part of the paper were obtained for the following 
parameters: i) The ground state was obtained by DMRG 
using a matrix size xdmrg = 400 and undergoing 10 
sweeps of two-site DMRG before switching to 40 runs 
for single-site DMRG. ii) The model consists of L = 150 
sites. Upon performing one of the three above described 
quenches (I, II or III) we used bias voltages in a range 
of Vb/A = (0,42). We always started from a half filled 
system in the canonical ensemble with total S z = 
and alternating up and down spins are chosen, iii) The 
time evolution was performed using a TEBD matrix 




FIG. 10: (Color online) Maximum simulation time reachable 
(QT I, xtebd = 2000) due to accumulation of entanglement 
entropy (left) and steady-state current (right). The left plot 
shows the time until a truncated weight of e c = 5 • 10 _J is 
reached at any bond of the chain for the first time. The right 
figure corresponds to the data in Fig. 2] Note the inverted 
color-scale in the left image, dark regions correspond to low 
values of the maximum time reachable. 



truncated weight of e c = 10~ 15 we dynamically adjusted 
the size of the TEBD matrices with a maximum matrix 
size of xtebd- We measured observables at each time 
step. 



Appendix B: Correlation of entropy and 
steady-state current 

The major limiting factor for time evolution using 
TEBD is the increase of bipartite entanglement^ 

Si = -tr (p L In (p L )) = -tr {p R In (p R )) , 

where Pl/ R denotes the reduced density matrix to the 
left (L) and to the right (R) of a bipartition at bond i. 
Using a maximum matrix dimension xtebd , we stop the 
simulation (for Fig. llO[) whenever the truncated weight 
at any bond exceeds a threshold value of e c = 5 • 10~ 5 , 
which defines our maximum simulation time (see 

Sec. HVAj) . In Fig-HUl wc pl Qt t e ] as a function of U 
and Vb (left) and compare it to the magnitude of the 
steady-state current for the same parameters (right). 
From our data we conclude that reachable simulation 
times due to accumulation of entanglement (and thus 
truncated weight) are non monotonic in U and Vb 
but can be characterized roughly by the magnitude 
of current in the system. We find this behavior to be 
generic to all investigated QTs. 
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where 



47vt,i I_YiL doJ ( (g £H+g£("))t' 2 -("- «/))(c.c)' 
Pl/r denote the electronic DOS of the (disconnected) left 
and right lead and g^/ H their advanced single particle 
Green's functions. 



